Trimmed and filtered reads were aligned using STAR (v2.7.0f). Read counting was performed using featureCounts (v1.6.2) with parameters “-s 1 -Q 255”. Gene annotation was generated by merging GENCODE (vM21), 3’UTR annotations in RefSeq (mm10, UCSC), and de novo 3’UTR sites discovered from the data.
Count table is read into R. A sample condition table is generated by parsing filenames.
# gene count file
counts_file <- "datasets/pit_utrseq_gencode_vM21_refSeq201906.txt"
# Obtain gene count table and sample conditions
temp <- get_count_table(counts_file, fixname = T)
count_table <- temp[[1]]
sampleCondition <- temp[[2]]
gene_info <- temp[[3]]
rm(temp)
colnames(count_table) <- gsub("p", "P", colnames(count_table))
# Get all chromosome Y genes and 3 chrX genes that are consistently female-biased
chrXYgenes <- as.character(gene_info[grepl("chrY", gene_info[,2]),1])
chrXYgenes <- c(chrXYgenes, "ENSMUSG00000086503.3", "ENSMUSG00000086503.4", "ENSMUSG00000099312.1")
chrMgenes <- as.character(row.names(gene_info[grepl("chrM", gene_info[,2]),]))
We used RUVseq to remove unwanted variables with empirically identified negative control genes.
# constructc edger expression set to obtain cpm values
y <- DGEList(counts=count_table, genes=row.names(count_table), group=sampleCondition$group)
y <- calcNormFactors(y)
count_cpm <- cpm(y)
# get spike-in gene names
spikes <- rownames(count_table)[grep("^ERCC", row.names(count_table))]
# filter for only genes expressed at cpm > 2 in at least 10 samples. Also removed mitochondrial genes and spike-ins
count_table_fil <- count_table[which(rowSums(count_cpm > 2) >= 10), ]
count_table_fil <- count_table_fil[!row.names(count_table_fil) %in% chrMgenes,]
count_table_fil <- count_table_fil[!row.names(count_table_fil) %in% spikes,]
# construct edger expression set again with filtered data
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$age)
y_fil <- calcNormFactors(y_fil)
# get log CPM counts
logCPM <- cpm(y_fil, log=TRUE, prior.count=1)
# Correct batch effect/unwanted variables using RUV seq
row.names(sampleCondition) <- colnames(count_table_fil)
set <- newSeqExpressionSet(as.matrix(count_table_fil),
phenoData = sampleCondition)
# upper quartile normalzation
set <- betweenLaneNormalization(set, which="upper")
# identify empirical control genes by performing an DE analysis across all conditions
design <- model.matrix(~group, data=pData(set))
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$group)
y_fil <- calcNormFactors(y_fil)
y_fil <- estimateDisp(y_fil, robust = TRUE, design = design)
fit <- glmQLFit(y_fil, robust = T, design)
lrt <- glmQLFTest(fit, coef=2:10)
top <- topTags(lrt, n=nrow(count_table_fil))$table
top$genename <- genename[top$genes]
# identified any genes with an FDR > 0.1 as negative controls
empirical <- subset(top, FDR > 0.1)$genes
# RUVseq normliazation
set1 <- RUVg(set, empirical, k=2)
# obtain normalized counts
logCPMc <- log2(normCounts(set1)+1)
qPCR of 183 puberty-related and control genes have been performed previously using the same samples. Here, we compared the expression quantification using 3’UTRseq (log2 normalized read counts) and that using qPCR (delta Ct)
# get qPRC results
all_data_normed_info_noC <- readRDS("datasets/fluidigm_all_data_normed_info_noC.rds")
# fix gene names
all_data_normed_info_noC$Primer[all_data_normed_info_noC$Primer == "Rab7l1"] <- "Rab29"
fluidigm_genes <- unique(all_data_normed_info_noC$Primer)
all_genes <- fluidigm_genes
# reformat pcr results table
PCR_puberty <- subset(all_data_normed_info_noC, tissue == "P") %>%
mutate(Sample = gsub("_", "", Sample)) %>%
group_by(Sample, Primer, age, gender) %>%
summarise(Ct= mean(Ct)) %>%
as.data.frame()
# compare pcr results and utrseq results (before and after normalization, only post-normalization data shown)
exp_compare_cor_values <- compare_qPCR(logCPM, logCPMc, PCR_puberty, outname, method = "spearman", ifgenplot = F)
compare_sum <- compare_qPCR_summaryplot(exp_compare_cor_values$cor_values)
exp_compare_cor_values$plots_after_correction$Pd37M4
Figure 1: Figure S1C: summary of UTRseq vs. qPCR quantification comparison
compare_sum[[2]]
Figure 2: Figure 1C: example scatter plot comparing UTRseq and qPCR
Pearson correlation coefficient between normalized UTRseq counts and qPCR delta Ct ranges from 0.7 to 0.78.
PCA plot and correlation heatmap using RUVseq normalized counts (log2 transformed).
# get PCA plots
logCPMc_PCAplot <- gen_pca_plot(logCPMc, npc=20, sample_label=F, label =F)
logCPMc_PCAplot[[1]]
Figure 3: Figure 1D: PCA plot of all samples
hmcol <- colorRampPalette(c("orange", "white", "darkmagenta"))(250)
plotHM(logCPMc, cols = hmcol, new_col_names = sapply(strsplit(colnames(logCPMc), "_"), "[[", 2), show_col_names = F, clust_method = "complete", anno = sampleCondition)
Figure 4: Figure S2A: correlation heatmap of samples
# calculate mean person correlation between biological reps
sample_cors <- cor(logCPMc)
sample_cons <- unique(sampleCondition$group)
x <- sample_cons[1]
sample_cors_means <- lapply(sample_cons, function(x) {
cors <- sample_cors[grepl(x, row.names(sample_cors)), grepl(x, colnames(sample_cors))]
cors[cors==1] <- NA
return(list(mean(cors, na.rm=T), cors))}
)
mean_cor <- round(sapply(sample_cors_means,"[[", 1),2)
Average pearson correlation between biological replicates (samples from the same biological condition) ranges from 0.95 to 0.97
Pairwise differential analysis was performed using EdgeR. Genes with a FDR < 0.05 and an absolute fold change> 1.5 are considered significant.
Comparisions are performed between male and female samples at each age
all_detected_genes <- genename[row.names(count_table_fil)]
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$group)
y_fil <- calcNormFactors(y_fil)
design <- model.matrix(~0+group+W_1+W_2, pData(set1))
y_fil <- estimateDisp(y_fil, robust = TRUE, design = design)
fit <- glmQLFit(y_fil, design=design, robust = TRUE)
# compare between sexes at each age.
pair_contrasts <- makeContrasts(d12_sex = groupd12F - groupd12M,
d22_sex = groupd22F - groupd22M,
d27_sex = groupd27F - groupd27M,
d32_sex = groupd32F - groupd32M,
d37_sex = groupd37F - groupd37M,
levels = design)
get_test_result <- function(contrasts=pair_contrasts, name){
qlf <- glmQLFTest(fit, contrast = contrasts[, name])
topTg <- topTags(qlf, n=nrow(y_fil$counts))
de <- as.data.frame(topTg[[1]])
de$genename <- genename[de$genes]
return(de)
}
de_result_list <- lapply(colnames(pair_contrasts), function(x) get_test_result(pair_contrasts, x))
names(de_result_list) <- colnames(pair_contrasts)
de_sig_list <- lapply(de_result_list, function(x) subset(x, FDR < 0.05 & abs(logFC) > log2(1.5)))
names(de_sig_list) <- colnames(pair_contrasts)
templist <- flatten(lapply(de_sig_list, function(x) list(subset(x, logFC>0), subset(x, logFC<0))))
names(templist) <- flatten(lapply(names(de_sig_list), function(x){
sapply(c("up", "down"), function(y) paste(x, y, sep = "_"))
}))
sig_genes_names_list_direction <- lapply(templist, "[[", "genename")
sig_genes_names_list <- lapply(de_sig_list, "[[", "genename")
sex_sig_genes_df <- get_sig_genes_df_from_list(de_sig_list)
colnames(sex_sig_genes_df) <- gsub("_sex_up", "_F", colnames(sex_sig_genes_df))
colnames(sex_sig_genes_df) <- gsub("_sex_down", "_M", colnames(sex_sig_genes_df))
upset(sex_sig_genes_df, sets=rev(c("d12_F", "d22_F", "d27_F", "d32_F","d37_F", "d12_M", "d22_M", "d27_M","d32_M", "d37_M")), keep.order = T,
order.by = "freq",
queries = list(list(query = intersects, params = list("d27_F", "d32_F", "d37_F"), color="tomato", active=T),
list(query = intersects, params = list("d27_M","d32_M", "d37_M"), color="steelblue", active=T)))
Figure 5: Figure 4A: overlap between sex-biased genes
# get sex biased genes in 2 ages
temp <- table(unlist(sig_genes_names_list_direction[c("d27_sex_up", "d32_sex_up","d37_sex_up")]))
f_biased_genes_used <- names(temp[temp !=1])
f_biased_genes_3 <- names(temp[temp ==3])
temp <- table(unlist(sig_genes_names_list_direction[c("d27_sex_down", "d32_sex_down","d37_sex_down")]))
m_biased_genes_used <- names(temp[temp !=1])
m_biased_genes_3 <- names(temp[temp ==3])
# f_biased_genes_used and m_biased_genes_used are used as input for Lisa
sig_df <- melt(sig_genes_names_list) %>%
separate(L1,into = c("age",NA))
# plot selected d12-d22 sex biased genes
ave_expr <- rowMeans(logCPMc)
names(ave_expr) <- genename[names(ave_expr)]
# figure2 genes are selected based on 1) d12 sex biased genes , 2) gene function mentioned in the text and 3) expression levels
figure2_genes1 <- c("Chrna4","Kcna4","Th","Drd4", "Asb4", "Fshb","Lhb","Serpine2", "Nupr1","Gpr101")
factor_plot(logCPMc, figure2_genes1, ncols = 5, sig_df = sig_df)
# get top m-/f-biased genes that are consistent from PD27
top_fbiased <- sort(ave_expr[names(ave_expr) %in% f_biased_genes_3], decreasing = T)[1:5]
top_mbiased <- sort(ave_expr[names(ave_expr) %in% m_biased_genes_3], decreasing = T)[1:5]
factor_plot(logCPMc, c(names(top_fbiased), names(top_mbiased)), ncols = 5, sig_df = sig_df)
enrich_list_direction <- lapply(sig_genes_names_list_direction, function(x) gprofiler(x, organism = "mmusculus", custom_bg = all_detected_genes, hier_filtering = "moderate", src_filter = c("GO", "KEGG","REAC","CORUM","HP","HPA","OMIM"),min_set_size = 5, max_set_size = 5000,min_isect_size = 3))
male_biased_pathways <- enrich_list_direction[c("d27_sex_down", "d32_sex_down", "d37_sex_down")]
female_biased_pathways <- enrich_list_direction[c("d27_sex_up", "d32_sex_up", "d37_sex_up")]
maleenrichmentmap <- enrichmap_like_pie(male_biased_pathways, layout_used = "kk", radius_scale = 0.05)
femaleenrichmentmap <- enrichmap_like_pie(female_biased_pathways, layout_used = "kk", radius_scale = 0.07)
maleenrichmentmap[[2]]
Figure 6: Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37
maleenrichmentmap[[3]]
Figure 7: Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37
femaleenrichmentmap[[2]]
Figure 8: Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37
femaleenrichmentmap[[3]]
Figure 9: Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37
# sex age interaction
sex_age_contrasts <- makeContrasts(d12_d22 = (groupd22F - groupd12F) - (groupd22M - groupd12M),
d22_d27 = (groupd27F - groupd22F) - (groupd27M - groupd22M),
d27_d32 = (groupd32F - groupd27F) - (groupd32M - groupd27M),
d37_d32 = (groupd37F - groupd32F) - (groupd37M - groupd32M),
d32_d22 = (groupd32F - groupd22F) - (groupd32M - groupd22M),
d37_d27 = (groupd37F - groupd27F) - (groupd37M - groupd27M),
d37_d22 = (groupd37F - groupd22F) - (groupd37M - groupd22M),
levels = design)
de_sex_age_result_list <- lapply(colnames(sex_age_contrasts), function(x) get_test_result(sex_age_contrasts, x))
names(de_sex_age_result_list) <- colnames(sex_age_contrasts)
de_sex_age_sig_list <- lapply(de_sex_age_result_list, function(x) subset(x, FDR < 0.05 & abs(logFC) > log2(1.5)))
sig_sex_age_genes_names_list <- lapply(de_sex_age_sig_list, function(x) x$genename)
#sig_genes_puberty <- lapply(sig_genes_names_list, function(x) x[x %in% all_puberty_genes])
d12_d22_sex_diff <- de_sex_age_sig_list$d12_d22$genename
#factor_plot(logCPMc, sort(d12_d22_sex_diff), ncols=5, name=F)
The genes shown in this heatmap includes:
that are not already shown in Figure 2
# get sex biased genes in 2/3 ages
sex_biased_genes_to_show <- unique(c(unlist(sig_genes_names_list[1:2]),d12_d22_sex_diff))
# exclude genes that are shown in figure 2 already
sex_biased_genes_to_show <- sex_biased_genes_to_show[!sex_biased_genes_to_show %in% c(figure2_genes1, names(top_fbiased), names(top_mbiased))]
#sex_biased_genes_to_show <- c(sex_chr_genes, "Kcna4","Th", d22_rest)
sex_genes_sum <- table(unlist(sig_genes_names_list[1:5]))
temp_df <- data.frame(value = sex_biased_genes_to_show, stringsAsFactors = F)
sex_gene_df <- left_join(temp_df, melt(sig_genes_names_list_direction)) %>%
# filter(value %in% sex_biased_genes_to_show) %>%
separate(L1, "_", into = c("age", NA, "dir")) %>%
spread(age, dir) %>%
mutate(n = sex_genes_sum[as.character(value)]) %>%
dplyr::select(value, d12, d22, d27, d32, d37) %>%
arrange(d12, d22, d27, d32, d37) %>%
column_to_rownames("value")
#sex_gene_df$d12_d22 <- ifelse(row.names(sex_gene_df) %in% d12_d22_sex_diff, "yes","no")
anno_cols_sex <-rep(list(c("up" = "tomato","down" = "steelblue")), 5)
names(anno_cols_sex) <- c("d12","d22","d27","d32","d37")
#anno_cols_inter <- list("d12d22" = c("yes" = "darkmagenta"))
p <- plot_gene_hm(row.names(sex_gene_df),
logCPMc,
clust_cols = F,
clust_rows = F,
row_cut = 1,
anno_colors = c(anno_cols, anno_cols_sex),
anno_rows = sex_gene_df[, c("d12","d22","d27","d32","d37")],
colnames = F,
# fontsize_row = 8,
colors = hmcol_yp,
sampleCon = sampleCondition)
plot(p$gtable)
Figure 10: Figure S3A: heatmap of selected sex-biased genes
These genes are TFs that are sex-biased at 2 out of 3 ages (PD27, 32, and 37) and are not shown in Figure 3.
# plots sex-biased tfs
sex_tfs <- c(m_biased_genes_used[m_biased_genes_used %in% alltfs], f_biased_genes_used[f_biased_genes_used %in% alltfs])
temp_df <- data.frame(value = sort(sex_tfs), stringsAsFactors = F)
sex_gene_df <- left_join(temp_df, melt(sig_genes_names_list_direction)) %>%
separate(L1, "_", into = c("age", NA, "dir")) %>%
filter(!is.na(age)) %>%
# spread(age, dir) %>%
pivot_wider(names_from = "age", values_from = "dir") %>%
mutate(sex = ifelse(d27=="down"|d32=="down"|d37=="down","down","up")) %>%
group_by(sex) %>%
arrange(d27, d32, d37,.by_group = T) %>%
column_to_rownames("value")
anno_cols_sex <-rep(list(c("up" = "tomato","down" = "steelblue")), 5)
names(anno_cols_sex) <- c("d12","d22","d27","d32","d37")
p_hm_tf <- plot_gene_hm(row.names(sex_gene_df),
logCPMc,
clust_cols = F,
clust_rows = F,
row_cut = 1,
anno_colors = c(anno_cols, anno_cols_sex),
anno_rows = sex_gene_df[, c("d27","d32","d37")],
colnames = F,
# fontsize_row = 8,
colors = hmcol_yp,
sampleCon = sampleCondition)
plot(p_hm_tf$gtable)
Figure 11: Figure 3B: Expression of sex-biased TFs
run_hyper <- function(query, geneset, total){
q <- length(intersect(query, geneset)) #num overlap geneset
m <- length(intersect(geneset, total)) # num of geneset
n <- length(total[!total %in% geneset]) # background that are not geneset
k <- length(intersect(query, total)) # num query genes
fold <- round(q/(k*m/(m+n)),3)
p <- phyper(q, m, n, k, lower.tail = F)
return(c("num_ol"=q,"num_genes_in_geneset"=m,"num_background_noneGeneset"=n,"num_query"=k,"fold_enrichment"=fold,"pval"=p))
}
# read in all disease genes
disease_genes <- read.table("datasets/pituitary_puberty_genes_combined_2020-09-28.txt", header = T, sep = "\t", stringsAsFactors = F, as.is = T)
all_disease_genes <- disease_genes$MGI.symbol
all_sex_biased_genes <- unique(unlist(sig_genes_names_list_direction[1:10]))
hyper_result <- run_hyper(all_sex_biased_genes, all_disease_genes, all_detected_genes)
hyper_result[c("fold_enrichment","pval")]
## fold_enrichment pval
## 1.801000e+00 2.667187e-05
# run CEMtools with default parameters
cem <- cemitool(as.data.frame(logCPMc), merge_similar = T, network_type = "signed")
# get module assignment
modules <- module_genes(cem)
modules$genename <- genename[modules$genes]
cemi_genes <- modules$genes
# transform into a list format
module_list <- lapply(split(modules, modules$modules), "[[", "genes")
module_list_genename <- lapply(split(modules, modules$modules), "[[", "genename")
# get average scaled expression for each gene for plotting in the next section
ave <- ave_count_table(logCPMc)
# get top hub genes
hubs <- CEMiTool::get_hubs(cem, n=20)
hubs_genename <- lapply(hubs, function(x) genename[names(x)])
-Left:Heatmap showing the scaled normalized expression levels of each gene, clustered separately within modules.
-Middle: plot summarizing scaled normalized expression patterns of each gene in each module. Thick line represent the average profile of the module.
-Right: expression profile of the top hub genes of each module.
-Bottom: top 10 hub genes for each module
# scale normalized count data (log transformed) and reorder columns
plotdata <- t(scale(t(logCPMc)))
colnames(plotdata) <- sapply(strsplit(colnames(plotdata), "_"), "[[", 2)
plotdata <- plotdata[, c(colnames(plotdata)[grepl("M", colnames(plotdata))], colnames(plotdata)[grepl("F", colnames(plotdata))])]
get_clust_gene_order <- function(clust_num){
# function to cluster genes within each cluster
used_genes <- module_list[[clust_num]]
used_data <- plotdata[row.names(plotdata) %in% used_genes, ]
phm <- pheatmap(used_data,
clust_cols = F,
clustering_distance_rows = "correlation",
clustering_method = "ward.D2",
silent = T)
order <- used_genes[phm$tree_row$order]
return(order)
}
# reorder genes based on clustering results within each cluster
reorder <- unlist(lapply(1:9, get_clust_gene_order))
cluster_length <- sapply(module_list, length)
# column annotation data
col_anno <- pData(set1)
row.names(col_anno) <-sapply(strsplit(row.names(col_anno), "_"), "[[", 2)
# set cutoffs for extreame values
all_values <- as.vector(as.matrix(plotdata))
custome_breaks <- unique(c(min(all_values), seq(quantile(all_values, 0.005), 0, length = 125), seq(0, quantile(all_values, 0.995), length=125), max(all_values)))
# generate heatmap
phm_modulegenes <- pheatmap(plotdata[reorder, ],
cluster_rows = F,
cluster_cols = F,
show_rownames = F,
annotation_col = col_anno[, c("age", "sex")],
annotation_colors = anno_cols,
color = hmcol,
breaks = custome_breaks,
gaps_row = cumsum(cluster_length),
silent = T)
# get a summarized plot of module gene expression
ave <- ave_count_table(logCPMc)
module_assignment <- modules$modules
names(module_assignment) <- modules$genes
pclusters <- plot_clusters(module_assignment, ave, "", ncols=1) + theme(axis.text.x = element_blank(), axis.text.y = element_text(size=14), strip.text = element_text(size=14))
# plot top module hub genes
top_hub_genes <- c("Tshb", "GH", "Prl", "Fn1", "Dlk1", "Ank1", "Mki67", "Col23a1", "Fshb")
tophub_genes <- factor_plot(logCPMc, top_hub_genes, ncols = 1)
hub_df <- as.data.frame(sapply(hubs_genename, function(x) paste(x[1:10], collapse = ",")))
colnames(hub_df) <- "Top10_hub_genes"
# put together
p <- as_ggplot(phm_modulegenes$gtable) + pclusters + tophub_genes +
plot_layout(widths = c(4,1.5,1)) +
theme_bw()
p + gridExtra::tableGrob(hub_df)
Figure 12: Figure 4:co-expression module expression and hub genes
top5_hub_genes <- unname(unlist(lapply(hubs_genename, function(x) x[1:5])))
factor_plot(logCPMc, top5_hub_genes, ncols = 5)
Figure 13: Expression patterns of top 5 hub genes
module_ME <- mod_summary(cem, "eigengene")
module_ME_df <- melt(module_ME) %>%
separate(variable, "_", into=c("id", "condition")) %>%
mutate(age=substr(condition, 2,4), sex=substr(condition, 5,5), rep=substr(condition, 6,6))
row.names(module_ME) <- module_ME$modules
module_ME <- module_ME[, -1]
p <- pheatmap(cor(t(module_ME)), silent = T)
plot(p$gtable)
Figure 14: Figure S5A: correlation of module eigengenes
ggplot(module_ME_df) +
geom_rect(xmin = -Inf,xmax = Inf,
ymin = -Inf,ymax = Inf,fill="gray90", alpha = 0.1) +
geom_bar(aes(x=id, y=value, fill=sex), stat="identity") +
facet_grid(modules~age, scales="free") +
# theme(axis.text.x = element_text(angle=45, hjust=1)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank()) +
scale_fill_manual(values = c("tomato", "steelblue"))
Figure 15: Figure S5B: module eigengenes
module_gene_enrich <- lapply(module_list_genename, function(x) gprofiler(x, organism = "mmusculus", custom_bg = all_detected_genes, hier_filtering = "moderate", min_set_size = 5, max_set_size = 5000, min_isect_size = 3, src_filter = c("GO")))
p1 <- plot_cluster_enrichment(enrich_list = module_gene_enrich[1:4], "CEMi_modules_enrichment_", width=13, exclude_domain = "hp", save = F)
p2 <- plot_cluster_enrichment(enrich_list = module_gene_enrich[5:9], "CEMi_modules_enrichment_", width=13, exclude_domain = "hp", save = F)
p1+p2 +
plot_layout(guides = "collect")
Figure 16: Figure S5C: pathway enrichment result of co-expression modules
save.image("mouse_pituitary_mRNA_analysis.RData")
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /opt/R/3.6.0/lib/R/lib/libRblas.so
## LAPACK: /opt/R/3.6.0/lib/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_CA.UTF-8 LC_COLLATE=en_CA.UTF-8
## [5] LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8
## [7] LC_PAPER=en_CA.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] pheatmap_1.0.12 CEMiTool_1.8.2
## [3] UpSetR_1.4.0 patchwork_1.0.0
## [5] scatterpie_0.1.4 ggraph_1.0.2
## [7] tidygraph_1.1.2 gProfileR_0.6.7
## [9] fgsea_1.10.0 Rcpp_1.0.1
## [11] dendextend_1.12.0 ggdendro_0.1-20
## [13] cowplot_0.9.4 gridExtra_2.3
## [15] reshape2_1.4.3 RUVSeq_1.18.0
## [17] EDASeq_2.18.0 ShortRead_1.42.0
## [19] GenomicAlignments_1.20.1 SummarizedExperiment_1.14.0
## [21] DelayedArray_0.10.0 matrixStats_0.54.0
## [23] Rsamtools_2.0.0 Biostrings_2.52.0
## [25] XVector_0.24.0 BiocParallel_1.18.0
## [27] ggrepel_0.8.1 ggpubr_0.2.1
## [29] magrittr_1.5 forcats_0.4.0
## [31] stringr_1.4.0 dplyr_0.8.3
## [33] purrr_0.3.3 readr_1.3.1
## [35] tidyr_1.0.0 tibble_2.1.3
## [37] ggplot2_3.2.1 tidyverse_1.2.1
## [39] pcaMethods_1.76.0 Biobase_2.44.0
## [41] gplots_3.0.1.1 RColorBrewer_1.1-2
## [43] rtracklayer_1.44.4 GenomicRanges_1.36.1
## [45] GenomeInfoDb_1.20.0 IRanges_2.18.1
## [47] S4Vectors_0.22.0 BiocGenerics_0.30.0
## [49] edgeR_3.26.5 limma_3.40.2
##
## loaded via a namespace (and not attached):
## [1] ggthemes_4.2.0 R.methodsS3_1.7.1 coda_0.19-2
## [4] ggpmisc_0.3.1 acepack_1.4.1 ffbase_0.12.7
## [7] bit64_0.9-7 knitr_1.33 aroma.light_3.14.0
## [10] R.utils_2.9.0 data.table_1.12.2 rpart_4.1-15
## [13] hwriter_1.3.2 RCurl_1.95-4.12 doParallel_1.0.14
## [16] generics_0.0.2 GenomicFeatures_1.36.4 preprocessCore_1.46.0
## [19] RSQLite_2.1.2 europepmc_0.3 bit_1.1-14
## [22] enrichplot_1.4.0 xml2_1.2.2 lubridate_1.7.4
## [25] assertthat_0.2.1 viridis_0.5.1 xfun_0.29
## [28] hms_0.5.2 jquerylib_0.1.4 evaluate_0.14
## [31] DEoptimR_1.0-8 fansi_0.4.0 progress_1.2.2
## [34] caTools_1.17.1.2 readxl_1.3.1 igraph_1.2.4.1
## [37] DBI_1.0.0 geneplotter_1.62.0 htmlwidgets_1.3
## [40] ellipsis_0.2.0 backports_1.1.4 bookdown_0.11
## [43] annotate_1.62.0 biomaRt_2.40.0 vctrs_0.2.0
## [46] withr_2.1.2 ggforce_0.2.2 triebeard_0.3.0
## [49] robustbase_0.93-5 checkmate_1.9.3 sna_2.4
## [52] prettyunits_1.0.2 cluster_2.0.8 DOSE_3.10.1
## [55] lazyeval_0.2.2 crayon_1.3.4 genefilter_1.66.0
## [58] pkgconfig_2.0.2 labeling_0.3 tweenr_1.0.1
## [61] nlme_3.1-139 nnet_7.3-12 rlang_0.4.2
## [64] lifecycle_0.1.0 modelr_0.1.5 cellranger_1.1.0
## [67] polyclip_1.10-0 graph_1.62.0 Matrix_1.2-17
## [70] urltools_1.7.3 base64enc_0.1-3 ggridges_0.5.1
## [73] viridisLite_0.3.0 bitops_1.0-6 R.oo_1.22.0
## [76] KernSmooth_2.23-15 blob_1.2.0 qvalue_2.16.0
## [79] robust_0.4-18 gridGraphics_0.4-1 ggsignif_0.5.0
## [82] scales_1.0.0 memoise_1.1.0 plyr_1.8.4
## [85] gdata_2.18.0 zlibbioc_1.30.0 compiler_3.6.0
## [88] intergraph_2.0-2 rrcov_1.4-7 cli_2.0.0
## [91] htmlTable_1.13.1 Formula_1.2-3 MASS_7.3-51.4
## [94] WGCNA_1.68 tidyselect_0.2.5 stringi_1.4.3
## [97] highr_0.8 yaml_2.2.0 GOSemSim_2.10.0
## [100] locfit_1.5-9.1 latticeExtra_0.6-28 GeneOverlap_1.20.0
## [103] grid_3.6.0 fastmatch_1.1-0 tools_3.6.0
## [106] rstudioapi_0.10 foreach_1.4.4 foreign_0.8-71
## [109] farver_1.1.0 digest_0.6.19 rvcheck_0.1.3
## [112] pracma_2.2.5 ff_2.2-14 broom_0.5.3
## [115] gRbase_1.8-3.4 httr_1.4.1 AnnotationDbi_1.46.0
## [118] colorspace_1.4-1 rvest_0.3.5 XML_3.98-1.20
## [121] splines_3.6.0 RBGL_1.60.0 ggplotify_0.0.3
## [124] fit.models_0.5-14 xtable_1.8-4 jsonlite_1.6
## [127] dynamicTreeCut_1.63-1 zeallot_0.1.0 R6_2.4.0
## [130] Hmisc_4.2-0 pillar_1.4.3 htmltools_0.3.6
## [133] glue_1.3.1 clusterProfiler_3.12.0 DT_0.7
## [136] DESeq_1.36.0 codetools_0.2-16 pcaPP_1.9-73
## [139] mvtnorm_1.0-11 lattice_0.20-38 network_1.15
## [142] gtools_3.8.1 GO.db_3.8.2 survival_2.44-1.1
## [145] rmarkdown_2.10.6 statnet.common_4.3.0 munsell_0.5.0
## [148] DO.db_2.9 fastcluster_1.1.25 GenomeInfoDbData_1.2.1
## [151] iterators_1.0.10 impute_1.58.0 haven_2.2.0
## [154] gtable_0.3.0